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Abstract 

We develop the strong-field approximation for high-order harmonic generation in hydrogen 
molecules, including the vibrational motion and the laser- induced coupling of the lowest two Born- 
Oppenheimer states in the molecular ion that is created by the initial ionization of the molecule. 
We show that the field dressing becomes important at long laser wavelengths (~ 2/im), leading 
to an overall reduction of harmonic generation and modifying the ratio of harmonic signals from 
different isotopes. 

PACS numbers: 33.80.Rv, 42.65.Ky 
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I. INTRODUCTION 



In the high-order harmonic generation (HG) process, an atomic or molecular system 
irradiated by intense laser light emits high-frequency coherent radiation. The properties of 
the emitted radiation have led to interesting applications in the past decade. To enumerate 



just a few, we mention the generation of coherent ultraviolet attosecond pulses 1 
the measurement of vibrational motion in molecules 4|, and tomographic reconstruction of 
molecular orbitals . This list shows that the main focus of HG has recently extended 

from atomic systems (mostly rare-gas atoms) to molecules, which possess more degrees of 
freedom, namely vibration and rotation. They enrich the dynamics and provide additional 
'control knobs' for HG. By understanding the effects of vibration and rotation on the HG 
spectra, one is able to manipulate the harmonic radiation. 

One of the main theoretical tools in understanding HG is the strong-field approximation 

n 

(SEA), also known as the Lewenstein model. Originally proposed to study HG in atoms [8|, 
it was later extended to molecules. It is the quantum-mechanical formulation of the three- 
step model 9j, which ascribes HG to a sequence of (i) ionization, (ii) acceleration of the 
continuum electron, and (iii) recombination. The three-step model has had great success in 
describing qualitatively the dynamics of harmonic generation. It also predicts correctly the 
value of the cutoff energy of the emitted harmonic radiation. Regarding the quantitative 
predictive power of the Lewenstein model for HG in molecules, we note that the model 
ignores the Coulomb forces acting on the active electron in the continuum. This affects 
most significantly the region of low harmonics, which is thus not accurately described. For 
high-order harmonics, the absolute value of the harmonic intensity is usually lower than the 
value obtained by numerically integrating the time- dependent Schrodinger equation (TDSE) 
(for atoms, see jlO| where such a comparison is made). Nevertheless, in the case of atoms, 
it was shown that the qualitative behavior of high-order harmonics usually agrees well with 
the TDSE result [l^. For molecules, the extra degrees of freedom can hinder the agreement 
with the exact results. One needs to consider different possible formulations and gauges to 
decide which one fits better the analysis of a given process. For an extended discussion about 
the choice of gauge and formulation in the context of molecules, see jll[ [l2]. For example, 
to describe the two-center interference effects [l3], it is advantageous to use the momentum 
formulation for the recombination step [12]. Based on our previous results [l2], we choose for 
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this work the length-gauge molecular Hamiltonian and the momentum formulation. More 
systematic, detailed quantitative comparisons to TDSE results are waiting to be performed. 

The vibration of the molecular ion, treated in the framework of the Born-Oppenheimer 
(BO) approximation, was included in 1^, 1^. Only one BO potential surface was taken 
into account, since at the wavelength of the commonly used Ti: Sapphire laser (800 nm), the 
other potential surfaces are expected to be irrelevant. The reason is that the laser field does 
not efficiently couple the BO surfaces of the molecular ion during the short time between 
ionization and recombination. At longer laser wavelengths, available for example at the 
Advanced Laser Light Source (ALLS) in Canada, the electron excursion times are longer 
and we expect that the inclusion of excited BO states is important. 

The physical picture we adopt is derived from the simple-man's model applied to 
molecules, the latter being directly linked to the physical interpretation emerging from 
the SFA. In this picture, after the active electron has reached the continuum, it does no 
longer interact with the core, until it returns to recombine. Meanwhile, vibrational motion 
takes place in the molecular ion. To describe this motion, we consider the two lowest BO 
potentials in the ion, coupled by the laser field. We investigate two different formulations: 
(i) the full SFA emerging from the integral equation for the evolution operator and (ii) a 
simplified model analogous to the atomic simple-man's model, but including the vibration 
and dressing of the ion. The simple semiclassical model (SM) succeeds in reproducing the 
general characteristics of the full SFA, while being much less demanding. 

In terms of computational effort, the newly proposed SFA model turns out to be very 
demanding. This is due to solving numerous times the time-dependent Schrodinger equation 
for the vibrational dynamics in the two coupled BO surfaces. The intense numerical effort 
restricts significantly the parameter space one can explore. To this end, we develop a saddle- 
point approximation, reducing the computational time dramatically, while accounting for 
the relevant physical mechanisms. The implementation of this technique is described in 
detail in the Appendix. 

The paper is organised as follows: in the first part of Section [Tll we introduce the new 
SFA model and briefiy discuss the computational details. The second part presents the 
extended SM model, giving an expression that quantifies the contribution of different electron 
trajectories (see pi|] and references therein) to the total HG spectrum. Section UTTl discusses 
the results, for three different laser wavelengths: short (800 nm), medium (1500 nm), and 
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long (2000 nm) wavelengths. The HG spectra for H2 and D2 as well as the ratios of harmonic 
intensities in the two isotopes are calculated. The characteristics of the results can be well 
understood from the analysis of the electronic trajectories in the framework of the SM 
model. The 2000 nm case fully shows the importance of taking into account more than 
one BO potential surface. Finally, in the last section we present our conclusions and future 
perspectives. Atomic units (a.u.) are used throughout this work, unless otherwise specified. 

II. THEORY 

A. The extended strong-field approximation model 

We consider high-order harmonic generation (HG) in the hydrogen molecule, allowing for 
molecular vibration to take place. The full Hamiltonian includes all Coulomb interactions 
between the electrons and the protons, but in our model the Coulomb repulsion between 
electrons is neglected except for using the correct ground-state energy and BO potential 
of H2. The direction of the molecular axis is kept fixed. In the presence of a linearly 
polarized electric field E(t), the length gauge molecule-field interaction is given by Hiit) = 
Hii{t) + Hi2{t), with Hiiit) = E(t) -ri and Hyiit) = Fj(t) ■r2. Here, ri and r2 are the electron 
coordinates. The integral equation for the full evolution operator U {t, 0) reads: 

U{t,0) = Uo{t,0)-t [ dt'U{t,t')[H-a{t') + H,2{t')]Uoit',0), (1) 

where Uq is the evolution operator for the field-free Hamiltonian. 

We calculate the total dipole acceleration a(t) as the time-derivative of the expecta- 
tion value Pdip of the dipole momentum pj| Pdip = —(Pi + P2)- Neglecting continuum- 
continuum contributions, it follows from Eq. ([T]): 

Pdip(t) ~ -^ f dt\^^''\t)\P^,^U{t,t')[Ha{t') + H,2{t'm'^°\t')) + c.c. (2) 
Jo 

Here, $™°^(t) = exp{—iEot) and Eq are the molecular ground state (including the 
vibrational coordinate) and its energy, respectively. Assuming that only one electron can be 
promoted in the continuum and neglecting the interaction of the continuum electron with 
the remaining ion, Eq. ([2]) simplifies to: 

Pdip(t) ^ 2z rrft'($-°i(t)|Pif/iv(t,t')?^2B(t,t')^ii(0l'^'"°'(i')) + C.C., (3) 

^0 
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with the evolution operators Uiy and IJ2B explained below. We used also the fact that the two 
electrons are equivalent and we neglected the exchange terms {3, [l8|. Within the linear 
combination of atomic orbitals (LCAO) approximation used in this work for describing the 
electronic orbitals of the initial molecule and of the molecular ion, the exchange term gives 
rise to an ionization matrix element that is zero for the transition to the electronic ground 
state of the ion, and non-zero for the transition to the first excited state. The contribution of 
the latter is negligible compared to the contribution of the direct term (which corresponds to 
transitions to the electronic ground state of the molecular ion), since the ionization potential 
for ionization leading into the excited state is bigger than that for the ground state of the 
molecular ion. 

Since the exact propagator is unknown, one must resort to approximations. Consequently, 
in Eq. (jSj) the full evolution operator has been factorized in a product of two operators: the 
Volkov propagator Uiy that describes the active electron in the continuum, neglecting the 
influence of the binding potential and the electron-electron interaction, and the operator U2B 
that propagates the wave function of the remaining molecular ion. The Volkov propagator 
oa. be decomposed spectrally, aud the resulting integration over momenta is further smrpl- 
fied by using the saddle-point method [8l]. For the propagator U2B we consider only the two 
lowest-lying BO potential curves. The two potential curves of the molecular ion correspond 
to the symmetric Cg state with electronic wave function and to the ant i- symmetric a-^ 
state with electronic wave function ip]^'^. Including dressing means that the two states are 
coupled by the external electric field via dipole coupling. At the same time, vibrational wave 
packets evolve in each BO potential. Finally, the dipole momentum reads: 



dip(t) ^ 2t [ dt' exp{-tS{p,,t,t')) 
Jo 



2tt 



1 3/2 



POO 

/ dRXo{R) Cl,{R, Ps, t) U2B{t, t') Cion{R, Ps, t') X0{R) + C.C., 

Jo 

with Ps = — j^i dt" A{t") / {t — t') being the saddle-point momentum and e a regularization 
parameter of the order of unity (we use e = 1). The quasiclassical action S is defined 
as 5(p,t,t') = jl,dt"[p + A{t")f/2 + \EQ\{t-t'). Here, A{t) = - f^^dt' E{t'). The 
action represents the phase accumulated by the free electron moving under the influence of 
the external field only (relative to the ground state). The vibrational ground state of the 
molecule is denoted by Xoi^)^ with R being the internuclear distance. 
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The ionization and the recombination steps are described by the transition matrix ele- 
ments Cion and £rec (which are vectors of c- numbers): 

£ion(i?,P,t) = (5) 

/ / dn dr2^{.(p, ri, r2)E(t) ■ riV'o-°i(i?, n, r2) 

//rfrirfr2V'o"°H^,ri,r2)Pi^v(p,ri,t)^r(^'r2) 

£rec(i?,P,t) = . . (6) 

/ / rfn rfr2 iPr iR^ ri, r2)Pi^v(p, ri, ^2) 

In Eqs. (E]) and ([6]), '?/'v(P5 ^) = exp(z(p + A(t)) ■ r)/(27r)^/^ is the spatial part of a Volkov 
solution with canonical momentum p of the electron, and ip^"^ is the electronic BO ground 
state of H2. The physical interpretation of the matrix elements in Eqs. ([5]) and (Ej) can be 
given in simple terms. They quantify how much of the initial vibrational wave function xo 
is transferred on each of the two BO potential surfaces of the molecular ion. The transition 
is made from the electronic ground state of the molecule to the intermediate state in the 
HG process via the dipole operator of the active electron. The intermediate state is the 
state in which the continuum electron is described by a Volkov state and the bound electron 
is in the or (Xg state of the molecular ion. The recombination matrix elements describe 
the recombination process of the active electron in the momentum form. We approximate 
the electronic ground state by ip^°^{R, ri, r2) = ip^°^{R, ri)ip^°'^{R, r2) and we use the LCAO 
approximation for and ip]^'^ (only Is hydrogenic functions are used, see Appendix). In 
this case, the lower matrix elements in Eqs. ([5]) and (EI) are identically zero, due to symmetry 
reasons. This would not be the case if another approximation, such as the Heitler-London 
wave function, was used. 

The evolution operator ?72B propagates the vibrational wave packets created by ioniza- 
tion in the two potential surfaces of the molecular ion, according to the time-dependent 
Schrodinger equation: 

f xTiR) \^( + yr(R) E(t) . D(i?) \ f xr(^) , 

where D(i?) is the transition dipole moment between the gerade and ungerade electronic 
states in the ion, Vg°^{R) are the ionic BO energy surfaces, and 171^ is the mass of one 
nucleus. The transition dipole moment points along the molecular axis. Its modulus can be 
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well approximated [19| by 

|D(i?)| = 0.4e-^ + i?/2. (8) 



The numerical propagation for Eq. ([7j) is described in detail in 20j. At the ionization 
time t' the initial vibrational wave packets x^g^{R) ^md Xu'^(-R) ^"^^ given by the initial 
vibrational ground state xo multiplied by the ionization matrix elements from Eq. ([5]). The 
propagation is carried out between the ionization time t' and the recombination time t in 
Eq. ([3]). Hereafter, we refer to the calculation that takes the dipole coupling between the 
energy surfaces into account as the two-level (2L) calculation. The case when the coupling 
is neglected is referred to as the one-level calculation (IL), analysed in detail in |21|. 

For large laser wavelengths, the calculation based on Eq. (jl]) becomes very time consum- 
ing. One solution is to employ the saddle-point method to approximate the integral over the 
ionization time t' . It gives the possibility to study a large parameter space for the laser field, 
while remaining close to the full SEA results and dramatically reducing of the computational 
time. The details of the saddle-point method are given in the Appendix. 

B. Simple-man's model including vibration 

In order to estimate the importance of field dressing in the process of HG, we investigate 



the electronic trajectories [22|, |23|] in the spirit of the simple-man's model [9j. For each of 
the classical electronic trajectories, we define a weight that assesses its contribution to the 
total harmonic spectrum. The trajectories are assumed to start with zero initial velocity 
at an arbitrary ionization time ti. Thereafter the electric field of the laser pulse accelerates 
the electron and at the moment when it returns to the core, the return kinetic energy is 
calculated. The return kinetic energy plus the ionization potential equals the photon energy 
of the emitted harmonic radiation. Between ti and tr the vibrational wave packet created 
in the molecular ion evolves in the two coupled potential surfaces (see Fig. [1]). For the case 
when only one level is taken into account in the molecular ion, the model was described 



m 



2l| . There, it was shown that the prediction of the model for the ratio of harmonics 
in D2 and H2 compares very well to the full SEA result at 800 nm laser wavelength. The 
simplified model has the advantage that it is easy to implement and it requires very little 
computational time compared to the full calculation, while containing all the significant 
physics. 
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FIG. 1: (Color online) Schematic view of the harmonic generation process in H2. Shown are the 
ground-state BO potential of the H2 and the two lowest BO potentials of H^. The three-step 
process consists of: (1) ionization, (2) evolution of the remaining molecular core as prescribed 
by the field-dressed ionic cxg and o"u potential-surfaces while the active electron is driven by the 
external field, and (3) recombination into the molecular ground state. 



in Eq. 



To derive the relevant equations, one^relies on applying the saddle point approximation 
, as done for the IL case in 



W{ti, tr 



exp 



2l|. The trajectory weight is: 
2(2/pp- 



3 l^(ti) 



Xo(^),0 U/2B(tr,ti) 



dRcos{krRcos6/2) x 







(9) 



where 6 is the orientation angle between the molecular axis and the polarization direction 
of the electric field E{t), Eq is the energy of the molecular ground state, and fcr is the return 
velocity of the electron k^{t^,t-^ = — j^' dt' A{t') / {t^ — t;) + A{t^). Equation ^ includes 
essentially all molecular effects on each electron trajectory: the instantaneous ionization 
rate at time tj, the motion of the vibrational wave packets on the coupled BO surfaces in the 
ion, and the 'cos' interference-term. The interference term appears due to the presence of 
the two molecular sites whose contributions to the harmonic radiation interfere. The wave 
packet spreading was not included in Eq. ([9]), in order to obtain a weight incorporating only 



the molecular effects. 
As described in 



22 



'he 'vertical' ionization potential Jp is defined in the Appendix. 

for electrons that tunnel out during the same optical cycle, 



there are two types of trajectories: the short trajectories (that last less than about three 
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quarters of the optical period) and the long trajectories. As it will be shown in the following, 
the influence of the molecular vibration on various trajectories is very different due to the 
different trajectory durations. 



III. RESULTS 



As explained in the beginning, we expect that the consequences of vibration and field 
dressing become more relevant at longer laser wavelengths. To this end, we analyse three 
different cases. For 800 nm laser wavelength the effects of field dressing are shown to be 
negligible. In the intermediate regime, for a laser wavelength of 1500 nm, the field dressing 
starts becoming important. Finally, at 2000 nm, the dressing can no longer be ignored. 
In our calculations, the laser electric field has a trapezoidal envelope, with 4 optical cycles 
turn-on and turn-off, and 6 optical cycles of constant amplitude. The laser field is linearly 
polarized, and unless specified, the molecule is oriented parallel to the laser polarization 
direction {i.e., we consider aligned molecules). The laser intensity is 5 x 10^^ W/cm^, unless 
stated otherwise. All results in this section have been obtained using the saddle-point 
approximation (see Appendix), since for long wavelengths, the computational times for the 
fully numerical SFA are prohibitive. Except in the range of low harmonics, the saddle-point 
IL calculations were found in excellent agreement with the harmonic spectra calculated by 
full numerical integration with the field dressing neglected; for more details, see Appendix. 
The propagation of the vibrational wave packets has been carried out on a grid. The grid 
method is slightly more accurate than the eigenfunction decomposition method employed in 



2l| for the XL SFA. 



A. 800 nm laser field 

Figure [2] compares the XL and the 2L calculations. The reader should ignore the region of 
harmonic orders below ^ X5. In the saddle-point approximation, the electronic trajectories 
with small travel times are not accounted for, so that the low-order harmonics are not treated 
accurately (see comment at the end of the Appendix). 

As expected, the difference caused by the field dressing is very small. This can be un- 
derstood, since the typical travel time in the harmonic generation process is of the order of 
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FIG. 2: (Color online) Harmonic generation using 800 nm laser pulses. Left column: IL case, right 
column: 2L case. The upper row shows the harmonic intensities for the H2 molecule. The lower 
row shows the harmonic ratio D2/H2 (black, continuous curve), and the ratio predicted if one takes 
into account the short trajectory only (red, dashed curve) or the long trajectory only (blue, dotted 
curve) . 

one optical cycle, during which the transfer of the vibrational wave packet on the excited 
potential surface of the molecular ion is negligible. A clear difference between D2 and H2 
appears when the ratio of harmonic intensities is taken (see bottom row). The ratio exhibits 
a strong variation around harmonic order 50. Comparing to the harmonic spectrum, this is 
the region where an interference minimum 13| in the harmonic emission occurs. Due to the 
different masses of the isotopes causing different speeds of nuclear motion, the position of 
the interference minimum slightly changes from one isotope to the other. The cutoff of the 
electronic trajectories approximates well the quantum-mechanical cutoff, as expected. For 
the trajectories, we show only the pair of short and long trajectories that has the biggest 
contribution to the harmonic spectrum. It is evident that the short trajectory is the one 
that reproduces well the exact ratio (see Fig. O bottom row) [21^ . We note that the part of 
the harmonic spectrum with frequencies below the value of the 'vertical' ionization potential 
(see Appendix) is not accessible for the classical trajectory analysis. 

To deeper understand the role of different electron trajectories in the harmonic spectrum. 



10 
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FIG. 3: (Color online) Trajectory weights for 800 nm laser pulses (H2 only). Left panel: the 
trajectory weights for the IL calculation. Right panel: the trajectory weights for the 2L calculation. 
The red, dashed curves are used for the short trajectory and the blue, dotted curves are for the 
long trajectory. The black circles correspond to the longer trajectories. The trajectories shown by 
the dashed and the dotted lines are used to calculate the ratios depicted by the same curves in 
Fig. El 

Fig. [3] shows their weights [see Eq. Q]. As it can be seen, different pairs of trajectories 
contribute to different regions of the harmonic spectrum. Each pair consists of a short and 
a long trajectory. The pair with the shortest excursion time was used to calculate the ratios 
shown in Fig. [3l since this pair has the biggest contribution. This explains its success in 
reproducing satisfactorily the exact ratio. The shorter trajectories are not strongly affected 
by the laser coupling at this wavelength, while the pairs with longer travelling time can be 
strongly affected by the field coupling (compare the lowest pairs in the left and right panel 
of Fig. [3l). These pairs do not contribute significantly to the total spectrum. 

B. 1500 nm laser field 

Due to the longer duration of an optical cycle for 1500 nm wavelength, one expects to 
see stronger signatures of the field dressing in the harmonic spectra. This is indeed evident 
in the results shown in Fig. |H Comparing the ratio of the harmonic signal (bottom row 
of Fig. H]) for the IL calculation and the 2L calculation, the short trajectory is again able 
to reproduce quantitatively the full results. In the 2L case, the ratio exhibits a somewhat 
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FIG. 4: (Color online) Harmonic generation using 1500 nm laser pulses. Left column: IL case, 
right column: 2L case. The upper row shows the harmonic intensities for the H2 molecule. The 
lower row shows the harmonic ratio D2/H2 (black, continuous curve), and the ratio predicted if 
one takes into account the short trajectory only (red, dashed curve) or the long trajectory only 
(blue, dotted curve). 



smoother variation with the harmonic order. We understand this by analysing the behavior 
of the trajectory weights in Fig. [51 Here, there is a significant influence of the field dressing on 
the trajectories. Namely, when the field coupling is included, the short trajectory is the one 
least affected, while the long one and the remaining less-contributing pairs become strongly 
damped. In such a case, the quantum interference effects between trajectories 22, 1231] have 
a smaller impact on the HG spectrum, so that the ratio becomes smoother. Again, note the 
strong variation of the ratio below harmonic order 150, which is related to the position of 
the interference minimum |l3|, as discussed above. 



C. 2000 nm laser field 

Finally, we study the case when the laser wavelength is large enough to allow the field 
dressing effects to fully manifest themselves. We show results for two laser intensities, in 
order to assess the influence of the fleld strength on the coupling of the BO surfaces. 

For the laser intensity 1=2.5 x 10^^ W/cm^, Fig. El compares the HG spectra for the IL 



12 



100 200 300 400 100 200 300 400 

Harmonic order Harmonic order 

FIG. 5: (Color online) Trajectory weights for 1500 nm wavelength (H2 only). Left panel: the 
trajectory weights for the IL calculation. Right panel: the trajectory weights for the 2L calculation. 
The red, dashed curves are used for the short trajectory and the blue, dotted curves are used for 
the long trajectory. The black circles correspond to the longer trajectories. The trajectories shown 
by the dashed and the dotted lines are used to calculate the ratios depicted by the same curves in 
Fig.H 

and the 2L case (top row). The reduction in the harmonic intensity is clearly visible on the 
scale of the graphs. Thus, one concludes that the field coupling can not be neglected at long 
wavelengths. 

The HG spectrum for the 2L case shows less trajectories interference and hence a smoother 
shape of the spectrum envelope . This feature is explained by the simple-man's analysis 
(Fig. [7]), which shows that the contribution of the longer trajectories to the HG spectrum is 
strongly attenuated due to the larger time the electrons spends in the continuum. For the 
same reason, the SFA harmonic ratio for the 2L case (bottom right panel of Fig. [H] ) agrees 
with the SM ratio much better than for the IL case (bottom left panel of Fig. [6]). In the IL 
case (left panel of Fig. [7]), the long trajectories with duration more than one cycle contribute 
significantly to the spectrum, at least for photon energies below the ~ 400**^ harmonic order. 
At higher harmonic orders, there is only one trajectory pair that contributes to the HG 
spectrum. Consequently, the SM model succeeds in reproducing well the full SFA ratio 
(bottom left panel of Fig. [H]). 

Both the IL and the 2L HG spectra show a minimum around the 120*^^ harmonic order. 
As discussed above, this minimum is caused by the interference of the emitted harmonic 



13 




100 200 300 400 500 100 200 300 400 500 

Harmonic order Harmonic order 



FIG. 6: (Color online) Harmonic generation using 2000 nm laser pulse with 2.5 x 10^^ W/cm^ 
intensity. Left column: IL case, right column: 2L case. The upper row shows the harmonic 
intensities for the H2 molecule. The lower row shows the harmonic ratio D2/H2 (black, continuous 
curve), and the ratio predicted if one takes into account the short trajectory only (red, dashed 
curve) or the long trajectory only (blue, dotted curve). 

radiation from the two molecular sites in D2 [3]. Noticeably, the minimum appears clearly 
in both the full SFA and the SM harmonic ratios. The presence of the cosine interference 
term in the expression for the trajectory weight given by Eq. is essential for reproducing 
this minimum. 

The effects of field coupling become even more apparent at higher laser intensity. In the 
following, we employ the value 1=5 x 10^^ W/cm^, used also for the shorter laser wavelengths. 
From Fig. [HI by comparing the HG spectra for the IL and the 2L cases, the difference between 
the IL and the 2L results is noticeable. 

In the IL case, the ratio is not well reproduced by the short trajectory, except beyond the 
700*^ harmonic. The reason is analogous to the lower-intensity case, as becomes clear when 
one analyses the trajectories weights shown in the left panel of Fig. [91 As more trajectories 
contribute to harmonic orders below 700, quantum interference takes place. Consequently, 
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FIG. 7: (Color online) Trajectory "weights for 2000 nm "wavelength and 2.5 x 10^^ W/cm^ intensity 
(H2 only). Left panel: the trajectory "weights for the IL calculation. Right panel: the trajectory 
"weights for the 2L calculation. The red, dashed curves are used for the short trajectory and 
the blue, dotted curves are for the long trajectory. The black circles correspond to the longer 
trajectories. The trajectories sho"wn by the dashed and the dotted lines are used to calculate the 
ratios depicted by the same curves in Fig. [6l 

one single trajectory is not enough to describe accurately the HG spectrum. In contrast, for 
harmonics of order higher than 700, there is only one short trajectory that dominates the HG 
spectrum. This clearly explains the agreement between the full result and the predictions 
of the SM model in this region. We have checked that these findings remains valid when 
additional factors due to wave-packet spreading are taken into account. 

The situation for the 2L case is strikingly different from the IL case. Namely, the ratio 
is much smoother, with reduced signs of interference. Such behavior is caused by the fact 
that most of the trajectories that would normally contribute to the spectrum become less 
important. The classical analysis shows that this is indeed the case (see Fig. [9]). The SM 
analysis shows that even the short trajectories are affected strongly by the field coupling. As 
a general characteristic of large wavelengths, the harmonic spectra tend the be 'smoothed 
out' by the field dressing {i.e., less interferences). In the IL case, for an extended part 
of the spectrum, there are many trajectories that contribute with comparable weights, so 
that the trajectory interference in the harmonic signal is strong. This makes the SM model 
inapplicable to explain the ratio. When the field dressing is included, the long trajectories 
are damped, hence the interference is almost absent. This reveals an important feature of 
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FIG. 8: (Color online) Harmonic generation using 2000 nm laser pulses with 5 x 10^^ W/cm^ 
intensity. Left column: IL case, right column: 2L case. The upper row shows the harmonic 
intensities for the H2 molecule. The lower row shows the harmonic ratio D2/H2 (black, continuous 
curve), and the ratio predicted if one takes into account the short trajectory only (red, dashed 
curve) or the long trajectory only (blue, dotted curve). 

the harmonic ratio, namely a maximum around harmonic order 600, which is due to second 
order destructive two-center interference in H2. 

In the remaining part of this section we consider the case when the molecule is perpendic- 
ular to the laser polarization direction {6 = 90°). The laser intensity is the 5 x 10^"^ W/cm^ 
as in Fig. [HI For this orientation, both dressing and two-center interference [13|] are absent. 
This is confirmed by the flat HG spectrum (lower part of Fig. [TUl) . The top row of Fig. [TU] 
shows the harmonic ratio obtained from the SFA compared to the SM ratio (left panel), 
and the trajectory weight (right panel). The SM model approximates the harmonic ratio 
well for harmonic orders higher than 700 for the same reason as in the IL case for parallel 
alignment (see trajectory weights in the upper-right panel of Fig. [TOj) . 

An important observation is that the shortest electron trajectory pair does not have 
the highest weight for the lower harmonics (upper-right panel of Fig. [TU]) . We have found 
numerically that the weight of the shortest pair increases with decreasing wavelength and 
becomes dominant at wavelengths lower than ^ 1300 nm. This behavior is related directly 
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FIG. 9: (Color online) Trajectory weights for 2000 nm wavelength and 5 x 10^'* W/cm^ intensity. 
Left panel: the trajectory weights for the IL calculation. Right panel: the trajectory weights for 
the 2L calculation. The red, dashed curves are used for the short trajectory and the blue, dotted 
curves are for the long trajectory. The black circles correspond to the longer trajectories. The 
trajectories shown by the dashed and the dotted lines are used to calculate the ratios depicted by 
the same curves in Fig. [HI 



to the temporal shape of the vibrational autocorrelation function 15|. It has a maximum 
at the time of the vibrational period of the ion, giving increased weight to long trajectories 
with suitable duration. Our results show that this effect plays a role only when the dressing 
is absent. 



IV. CONCLUSIONS 



In this work, we have analysed the possibility to include field dressing in the strong- 
field approximation for harmonic generation in H2 molecules. Previously, the vibration of 
the molecular ion formed upon ionization was considered to take place on the lowest BO 
potential surface only [21j. Here, we take into account two potential surfaces, coupled by 
the external field via the dipole interaction. Such a modification proves to be essential at 
long laser wavelengths. 

To study the effect of field dressing, we use an extension of the Lewenstein model which in- 
cludes fully quantum-mechanical vibrational motion of the molecular ion in the field-coupled 
lowest two BO surfaces. The resulting expression for the electronic dipole momentum turns 
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FIG. 10: (Color online) Harmonic generation using 2000 nm laser with 5 x 10^^ W/cm^ interacting 
with molecules perpendicular to the laser polarization direction (6 = 90°). Top row: Harmonic 
ratio D2/H2 compared to the SM ratios for shortest trajectories (left) and the trajectory weights for 
H2 from the simple-man's model (right). The red, dashed curves are used for the short trajectory 
and the blue, dotted curves are for the long trajectory. The black circles in the top-right panel 
correspond to the longer trajectories. Bottom: harmonic spectrum for H2. 

out to become numerically demanding even at the moderate wavelength of 800 nm. The rea- 
son is that a large number of one-dimensional time-dependent Schrodinger equations have to 
be solved. At longer wavelengths, the calculation becomes seriously prohibitive. Therefore, 
we applied the saddle-point method successfully, replacing the time integration by a sum 
over a few relevant terms. 

We have investigated three different laser wavelengths, 800, 1500, and 2000 nm. For the 
800 nm laser field, the effects of the laser dressing can be safely neglected. At 1500 nm 
wavelength and more strongly at 2000 wavelength, the effects of the field dressing manifest 
themselves in the harmonic spectrum, and more prominently in the ratio of harmonic in- 
tensities for D2 and H2. The field coupling has the effect of 'smoothing' the interferences in 
the harmonic spectrum, by lowering dramatically the contribution of the long trajectories 
to the HG spectrum. Consequently, only the shortest trajectory contributes significantly to 
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the spectrum. At 2000 nm, we found that this effect leads to a much clearer observation of 
two-center interference in the ratio D2 vs. H2. 

The proposed model could prove its usefulness in interpreting the experimental data avail- 
able from the use of the recently available long laser-wavelength sources and for uncovering 
new properties of the harmonic radiation in molecular systems. 

This work was supported by the Deutsche Forschungsgemeinschaft. 



APPENDIX: SADDLE-POINT APPROXIMATION 



The saddle-point method for the study of harmonic generation has been successfully 
applied in the context of the strong-field approximation {e.g., for atoms see 22], and for 



non-vibrating molecules see 



24|). due to the strongly oscillating phase factor exp(— ■iS') in 



the expression of the dipole moment [see Eq. (jlj)]. In the present work, we apply the saddle- 
point method to approximate the integral over the ionization time t' in Eq. (jl]). Compared to 
atoms and non- vibrating molecules, one encounters the additional problem that the part of 
the integrand which describes the molecular vibration can also oscillate with the integration 
variable. This oscillation cannot be calculated analytically. If the oscillatory vibrational 
part was known analytically, one could calculate the saddle points exactly. Based on our 
previous analysis of the IL case |2l|], we concluded that it is possible to isolate the desired 
oscillatory part. In order to do this, we make use of the fact that the main contribution to the 
electronic dipole moment comes from transitions between the electronic ground state of the 
initial molecule and that of the molecular ion. Consequently, we subtract from the energy 
curves of the molecular ion the quantity 6V = V^°'^{R), with R the average internuclear 
distance in the vibrational ground state of the neutral molecule. This means that, the 
propagation in Eq. (JTj) is done with the energy curves re-defined as V^'^'^ = V^'^'^ — 6V, while 
a term {t — t')6V has to be added to the action in Eq. (jlj). Combined with the ground-state 
energy Eq of the molecule, an effective vertical ionization potential Ip = \Eq\ + 6V appears 
in the semiclassical action. With this slight modification, the main oscillatory behavior is 
now concentrated in the re-defined semiclassical action, and one can safely apply the saddle- 
point method. There is still one difficulty: the LCAO approximation used to approximate 
the electronic states in both the molecule and the ion gives rise to ionization matrix elements 
in Eq. (j5j) that have a pole close to the saddle-point, since the upper ionzation matrix element 
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from Eq. ([5]) is proportional to 



E(t) ■ R ^.^ f[p + A{t)]-R 



[p + A(t)]-R 



2E(t) ■ [p + A(t)]cos 



2 2 
2 2 



-2 



+ 
-3 



(A.l) 



While for a given time t at which the electronic dipole moment is to be calculated, the 
equation for the saddle point t'^ reads: 



0, 



(A.2) 



the equation for the pole t'^ in the ionization matrix element corresponding to the transition 
to the gerade state of the ion [see Eq. flA.ip ] is: 



[Ps(tp + A(t; 



/ \12 



0. 



(A.3) 



The parameter Z in Eq. flA.Sp is the nuclear charge of the hydrogenic orbitals used in the 
LCAO approximation. We use Z = 1. The corresponding value for Z^/2 = 0.5 in Eq. (1A.3I1 
is close to and smaller than the value of Jp = 0.59 in Eq. flA.2p . As a consequence, the 
two critical points are very close to each other, with the pole having a smaller imaginary 
part than the saddle point. In this case, the usual saddle-point formula has to be modified 
accordingly to take into account the presence of the pole. For clarity, let us calculate the 
contribution of one saddle point and one paired pole Xp to the integral: 

ft 



(A.4) 



such that S'{xs) = and f{x) has a pole of order 3 in Xp. (According to Eq. flA.ll) . one has 
also a second-order pole. Its treatment is entirely similar to that of the 3rd-order pole.) To 
this end, we re-write the integral as 



dx f{x) 



{x - Xp)3' 



(A.5) 



with f(x) = f{x){x — Xp)^ bounded at the pole. According to the method described in 



25l . |26|, the first-order asymptotic contribution of the pair of critical points is 



-iS{xs 



dx 



r I /s /p / _ 

Jp \ \ X Xr 



Xa Xr 
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2 V*^ '^S 



[X — Xr 



(A.6) 
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FIG. 11: Comparison between the full calculation (black curves) based on Eq. ^ and by using 
the saddle-point approximation (red curves). The calculations are for the IL case. 



The integrals appearing in Eq. flA.6|) are of the form 



dx 



exp(— x^/2) 



[x - Xo 



(A.7) 



where xq 



Xr 



Xa 



Im(xo) < 0, and k is an integer. Equation (1A.7I) is adapted from 



the case with Im(a;o) > appearing in [27|. In Eq. flA.Tp . Dn{x) is the parabolic cylinder 
function, for which highly accurate numerical routines are available 28[. Thus, to calculate 
the dipole momentum from Eq. (jl]) at a given time t, one needs to find all saddle points 
and poles with the real part smaller than t and positive imaginary part, and calculate their 
contribution according to the above. In addition, the contribution from the pole has to 
be added, according to the Cauchy's theorem. The residue in the pole Xp is most simply 
calculated using directly the simplified expression in Eq. (lA.6p . which holds in the vicinity 
of the saddle point (where the pole is situated). 

The accuracy of the approximation can be seen in Fig. [TTl for various wavelengths. The 
saddle-point method gives more accurate results, the higher the laser wavelength, which can 
be seen from the figure. The comparison between exact and saddle-point calculation was 
carried out for the IL case using the eigenvalue decomposition {21 ] for numerical propagation 
of the vibrational wave packets, since in this case the calculation is much faster. The 
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eigenvalue decomposition, which can be used only for the IL case, has been applied only for 
the purpose of this comparison, while in the rest of this work we use numerical grids. For one 
laser wavelength, we also checked that the agreement with the saddle-point approximation 
holds as well for the 2L case (800 nm, not shown here). The harmonic-energy region where 
the agreement should be sought does not include the low-frequency harmonics. The reason 
is that all the saddle points close to the end point t of the time integration-interval are 
ignored in our approach. This is because their contribution to the integral can no longer 
be described by the above procedure. On the other hand, the critical points that are close 
to the recombination time t will contribute only to the low-frequency part of the harmonic 
spectrum, which is not the focus of this work. In practice, we ignore the saddle points for 
which Re(t — t'^) < Tq/IO , with Tq the period of one optical cycle. With this in mind, the 
saddle point method gives very accurate results for the high-energy part of the harmonic 
spectrum, while reducing significantly the computational time {e.g., for a laser wavelength 
of 1500 nm, by a factor of 20). 
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